Clustering and non-gaussian behavior in granular matter 
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We investigate the properties of a model of granular matter consisting of N Brownian particles on 
a line subject to inelastic mutual collisions. This model displays a genuine thermodynamic limit 
for the mean values of the energy and the energy dissipation. When the typical relaxation time r 
associated with the Brownian process is small compared with the mean collision time r c the spatial 
density is nearly homogeneous and the velocity probability distribution is gaussian. In the opposite 
limit t t c one has strong spatial clustering, with a fractal distribution of particles, and the 
velocity probability distribution strongly deviates from the gaussian one. 
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In the past few years granular materials have become 
an intriguing subject of research fjj - 0, since they pose 
novel questions and challenges to the theorists and exper- 
imentalists. The constituting elements of such materials 
are solid particles, whose size may range from few mi- 
crons to few centimeters, and which are subject to non- 
conservative contact forces such as friction and cohesion. 

Their collective behavior is peculiar and different from 
other forms of matter, such as solids, liquids or gases, 
and the ordinary statistical mechanical approach, which 
successfully deals with large assemblies of microscopic 
particles is not adequate. 

Generally speaking granular materials cannot be de- 
scribed as equilibrium systems neither from the config- 
urational point of view nor from the dynamical point 
of view. It is known in fact that these systems remain 
easily trapped in some metastable configurations which 
can last for long time intervals unless they are shaken 
or perturbed 0. On the other hand while in equilib- 
rium statistical mechanics the kinetic energy per particle 
is proportional to the temperature and the velocities are 
gaussianly distributed, in the systems we consider the 
tails of the distribution deviate from the Maxwell law 
H . This phenomenon is accompanied by a pronounced 
clustering of the particles or inelastic collapse 

Several approaches have been proposed for the study 
of the so-called "granular gases" [^J|] . One crucial differ- 
ence between ordinary gases and granular media is rep- 
resented by the intrinsic inelasticity of the interactions 
among the grains, which makes any theory based on en- 
ergy conservation , e.g. for ideal gases, not suitable. 

In the present work we study a one dimensional me- 
chanical model, in the spirit of the one recently in- 
troduced by Kadanoff and coworkers ||, but contain- 
ing some important differences regarding the energy- 
exchange process. We consider N identical particles on a 
circle of length L fnl obeying to the following equations: 
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where, 1 < i < JV, Tp is the temperature of a mi- 
croscopic medium that we discuss below, t is the re- 
laxation time, in absence of collisions, and fi(t) is a 
standard white noise with zero average and variance 
< fi(t)fj(t') >— 5ij5(t — t'). In addition the particles 
are subject to inelastic collisions according to the rule 
v[ — v'j = —r(vi — Vj), where r is the restitution coeffi- 
cient (r = 1 for the completely elastic case) fll| . 

The introduction of the viscous term takes into ac- 
count important factors, generally disregarded in simpli- 
fied models, namely the friction among particles and en- 
ergy transfers among different degrees of freedom, which 
are relevants in real granular systems. The damping 
and noisy terms are very natural when interactions be- 
tween particles and the environment (particle-fluid inter- 
actions) start to be important. Another important class 
of phenomena in which a viscous damping and a noisy 
term are naturally present is represented by the fluidized 
beds, where the vibration of the bottom of the box pro- 
duces a random force on the particles Q , j^] and p2[ . 

The model above differs from the one proposed by 
Kadanoff and coworkers because in the latter the par- 
ticles, subject to inelastic collisions, are confined to an 
interval of length L by two asymmetric walls: the first 
reflecting them elastically and the second supplying en- 
ergy to the particles according to a gaussian distribution 
at fixed temperature. As evident from the simulation 
of ref. H one observes a somehow trivial clustering of 
the particles next to the elastic wall. We found numeri- 
cally a more serious shortcoming of such model consist- 
ing in the fact that the average energy and the aver- 
age energy dissipated per particle, defined respectively 
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as the time average of E(t) = l/2Xii = i ^(i) 2 /^ and 
W{t) = (E(t 2 ) - E(ti))/(t2 - h) where t 2 and t\ rep- 
resent the times at which two successive collisions take 
place, are not independent from the total number of par- 
ticles, but decay exponentially as ~ exp(— cN), showing 
that the system does not possess a proper thermody- 
namic limit. The existence of thermodynamic limit in 
real granular systems is not clear, but it seems a rather 
natural requirement in a statistical mechanical approach. 

Due to the inelastic collisions, in order to reach a sta- 
tistically stationary situation some energy must be in- 
jected into the system. This is achieved in our model 
by the random noise term acting on each particle. This 
term mimics the action of a vibrating box. Notice that 
the boundary conditions and the energy-pumping mech- 
anism are different from that of reference || and present 
the advantage of providing a 'good' thermodynamic lim- 
its as far as the energy E and energy dissipation W are 
concerned, i.e. < E > and < W > become independent 
of N for large values of N. 

Moreover our system does not have walls and the clus- 
tering is non-trivial. On the other hand, the energy feed- 
ing mechanism adopted in forces to introduce a some- 
how artificial cooling of the particles, by renormalizing 
the velocity of the center of mass at each time step. In 
our formulation this procedure is overcome by the pres- 
ence of the thermal bath. 

For each given choice of r and r the system, after a 
long transient, reaches a stationary state with certain 
properties. The presence of two time scales, namely r and 
the mean collision time r c , leads to different dynamical 
regimes. 

a) When r <C t c it is easy to argue that the grains reach 
a rather simple statistical equilibrium and that their ve- 
locity distribution is that of an ideal gas with an effective 
temperature T F , slightly lower than the temperature of 
the 'heat bath', T F . 

b) In the opposite limit r 3> r c the driving mechanism 
towards the macroscopic stationary state is dominated 
by the collision process itself. 

Two phenomena are observed: 

1) the velocity distribution ceases to be gaussian and 
the deviation becomes more and more pronounced with 
decreasing values of the restitution coefficient r. 

2) The spatial distribution becomes strongly inhomo- 
geneous. 

It is worth to stress that, at variance with the clus- 
terization in Ref. ]^], in our case the clusters are created 
and destroyed continuously in the system as long as the 
system evolves. The inhomogeneity in the spatial distri- 
bution of the grains, see Fig. ([!]), can be quantitatively 
characterized by the so called Grassberger-Procaccia di- 
mension d 2 that we compute from the correlation func- 
tion C{R) defined as: 



(3) 

where T represents the duration of the simulation. In 
Fig.(^) we show C(R) vs R for a clustering situation. It 
turns out that d 2 is lower than 1 when r < 1, e.g. for 
t = 100 and r = 0.6, d 2 = 0.59, while for r = 100 and 
r = 0.9 d 2 ~ 1. 

Notice that in this case the stationary regime is 
brought about by the collisions and that these occur more 
frequently in the regions of higher density. For instance 
for r = 100 and r — 0.7 one gets for the number of colli- 
sions as a function of the spatial density p, N co ll (p) ~ P 2 ■ 

Assuming that in the stationary state the power dis- 
sipated through the collisions must balance the power 
adsorbed from the heat bath, one derives the following 
relation between the temperature Tp, the average energy 
E, the power dissipated by the collisions W and r fl3|| : 

W=-(T F -2E) (4) 

T 

Notice that in the absence of collisions (or r <C r c ) Tp = 
2E and W = as in the ideal gas. Eq. (|j) is well satisfied 
for different values of r and r. 
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FIG. 1. Snapshot of the particle density at a given time 
for T F = 1, r = 0.6, r = 100, N = 200, L = 200. All 
quantities are in arbitrary units. The dashed line represents 
the homogeneous density. 

As we quoted already in the inelastic regime one ob- 
serves a strong deviation from the gaussian behavior for 
the velocity distribution. In Fig.(0) we display the ve- 
locity distribution in a nearly elastic case (r = 0.01 and 
r = .99) and in a strong inelastic regime (r = 100 and 
r = 0.7). As it is possible to see it exists an evident de- 
parture from the gaussian behavior in the inelastic case, 
where the velocity distribution shows almost exponential 
tails. The above result is close to that observed in refs. 
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FIG. 2. C(J?) against 7? for two different choices of param- 
eters: t = 100, r = 0.6 (top) and r = 100, r = 0.99 (bottom). 
In both cases T F = 1, N = 200 and L = 200. The dimension 
takes respectively the values fib = 0.59 and g?2 = 1. 

Let us now try to relate the clustering properties of the 
system to the velocity distribution. In order to do that 
we consider the following quantities: the distribution of 
boxes, Ni, ox (m), containing a given number, m, of par- 
ticles and the average kinetic energy, Eki n {m), in a box 
occupied by m particles . Making the hypothesis that 
in each box the average velocity of the particles is zero, 
i.e. < v(m) >— (very well confirmed by the numerical 
data) , one finds that Eki n {m) provides a measure of the 
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FIG. 3. Rescaled velocity distribution P(v/a) against v/a: 
crosses are simulation data with r = 100, r = 0.7, the dashed 
line represents toy-model data, X's are simulation data with 
r = 0.01 and r = 0.99 and the dot-dashed line represents the 
gaussian distribution. In both cases Tf = 1, N = 200 and 
L = 200. 

variance of the velocity distribution in each box: 
Ekin{m) — \ < v 2 (m) >= icr 2 (m). We consider first 
the non-clusterized case (r << 1 and r ~ 1). Within 
this regime we find from the simulations that: 

X m c rn 

a 2 elas (m) ~ const. N^ a x s (m) = — — , (5) 



where A = N/Ndiscr is the average number of particles in 
each box and N^ s (m) is a Poisson distribution. By as- 
suming in each box a gaussian velocity distribution with 
a constant variance cr 2 las (m) it turns out that the global 
velocity distribution P e i as (v) is gaussian. Let us recall 
that the Poisson distribution is the one associated with 
a process of putting independently XN particles into N 
boxes. 

Let us turn to the non-elastic case. If r = 100 and 
r = 0.7, considering the occupied boxes (m > 0), we 
obtain from the simulations: 

p — am 

aidM-m-" NlZ l {m) = , (6) 

m 

with (3 ~ 0.5 and a ~ 0.14. Let us compute from 
these scalings the global velocity distribution. Taking 
into account that the spatial probability distribution of 
the particles is N^ l (m) and assuming that their lo- 
cal velocity distribution is gaussian, but with a variance 
af nei (m) ~ m _/3 which depends on the occupancy, we ob- 
tain, for the global velocity distribution Pi ne i{v), which 
in the continuum limit should correspond to 

P inel (v) ~ J2 e(-^e- am . (7) 

m— 1 

We stress how it exists an astonishing agreement between 
the numerical results and the ones obtained with a toy 
model which just makes the following hypothesis: 

i) non-Poissonian distribution for the box occupancy; 

ii) gaussian distribution of velocities in each box with 
a density-dependent variance. 

The hypothesis about the scaling relation between the 
velocity variance and the local density, apart from be- 
ing justified numerically, can be understood in the fol- 
lowing way. The stationarity and the scale-invariance 
of the cluster distribution, implies a certain distribution 
of lifetimes for the clusters. In particular each cluster 
has a lifetime which is inversely proportional to its size. 
The scale-invariant cluster-size distribution thus implies 
a scale-invariant distribution for the lifetimes. The clus- 
ter lifetime (its stability) is strictly related to the vari- 
ance of the velocity distribution inside the cluster itself. 
In order to ensure the stability of a cluster in a stationary 
state we have to require that the velocities of the parti- 
cles belonging to it are not too different, or equivalently 
that the variance of the distribution is smaller the higher 
the density. So, given a scale-invariant distribution of 
clusters one would expect a scale-invariant distribution 
of variances (^). An independent check is provided by 
the behavior of the average relative velocities in boxes 
with different numbers of particles. For r = 100 and 
r = 0.7 one obtains: 

w£fV)~m-'* (8) 

with 7 ~ 0.3, which indicates that the stability of a clus- 
ter is connected with the smallness of the velocity fluctu- 
ations inside it. In other words, the spatial clusterization 
corresponds to a clusterization in velocity space. 
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This analysis shows that in the present model it does 
exist a relation between the clustering phenomenon (in 
physical space or in velocity space) and the non-gaussian 
distribution of the velocities. 

A natural question to ask is whether our findings can 
be an artifact of the one dimensional dynamics. Can 
one expect to observe the non-Maxwell distribution and 
the clusterization even in higher dimensions? In order 
to clarify such issue, we consider a different stochastic 
process. All the particles perform the Brownian motion 
described by eqM if te[KAt, (K + l)At], (where K is 
an integer number) while at the instants tx = KAt each 
particle may collide with probability p with one of those 
particles which are spatially close to it according to the 
Boltzmann Stosszahlansatz (BS). 

In pra ctice we perform the following algorithm a la 
Bird pq| : at each discrete time tx for each particle, i, 
we extract out of a uniform distribution in the interval 
[0, 1] a random number y. If y > p there is no collision, 
otherwise the particle i scatters with another particle j 
if \xi(tx) — Xj(tfc)\ < I, and their collision probability is 
proportional to \vi(tx) — Vj{tx)\- This process renders 
the BS approximation exact; in fact in the limit N — > oo, 
p — > 0, I —> 0, At — > one can write for the stochastic 
process described above the following Boltzmann equa- 
tion for the one particle distribution P(x, v, t) fl6f| : 
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where A ~ p/At ~ l/r c and depends on the mean 
particle density. 

It is easy to show that in the limit of elastic collisions 
(r = 1) the velocity distribution becomes gaussian for 
any value of r, At and p. On the contrary for r ^ 1 an 
analytical solution is not known. 

In order to further clarify the relevance of the dimen- 
sionality, we have performed some simulations in 2d. In 
this case we find non gaussian tails for the velocity distri- 
bution, together with clusterization, at large values of r 
and strong inelasticity. In 2d the velocities in the collision 
change according to the following rule: 

Vi(t K + 0+) - Vj(t K + 0+) = re(vi(t K ) - v 3 {t K )) (11) 

where e is a two-dimensional unit vector of random ori- 
entation. In this case the results are qualitatively similar 
to the ones obtained in one dimension, i.e. at large values 
of r and strong inelasticity the Grassberger-Procaccia di- 
mension is smaller than 2 and the velocity distribution is 
not gaussian. 

In summary we introduced a model of granular gas 
with inelastic collisions between the particles. The sys- 
tem exhibits a variety of regimes ranging from a com- 
pletely elastic case without clusterization and gaussian 



distribution for the velocities to an inelastic regime with 
strong clusterization and non-gaussian velocity distribu- 
tion. In this framework we have shown a possible sce- 
nario to relate the clustering properties of the system to 
the velocity distributions. These results seem promising 
and give us the hope to be able to perform analytic work 
based on the Boltzmann approach in order to clarify fur- 
ther the model. 
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